# -*- coding:utf-8 -*-
'''
2021/04/30 

This program is created by ChenghaoQian in Zhejiang University, Hangzhou, Zhejiang, China
Email: chenghao_q@zju.edu.cn  chenghao_q@qq.com


'''
import numpy as np
import math
import random
import copy
import os
cur_path = os.getcwd()
os.chdir( cur_path )
import pandas as pd
import time
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


# mpl.use('qt5agg')


'step 0: set global variables'
#switch to the operating folder
cur_path = os.getcwd()
os.chdir( cur_path )
        
#set global size
x_size = (110) *10**-3 #尺寸，单位m
y_size = (110) *10**-3
z_size = (50)  *10**-3
x_coord_start = -0.5 * x_size
y_coord_start = -0.5 * y_size
z_coord_start = 0
grid_dist = 1E-3 #每个正方体cell的边长，单位m

#set delta time
dt = 0.5E-3 #单位s

C0 = 2.0 #速度更新的常数
C1 = 2.0 #混合分数更新的常数

Ka_crit = 0.1#varible need to be modified    可能0.05比较好
Pi_crit = 1# %2

init_float = np.nan   
    
class Coord(): #坐标类，单位mm
    def __init__(self):
        self.x = init_float
        self.y = init_float
        self.z = init_float
    
class Velo(): #速度类，单位m/s
    def __init__(self):
        self.x = init_float
        self.y = init_float
        self.z = init_float

class Mixfrc(): #混合分数，0-1间的数
    def __init__(self):  
        self.value = float(-1)
        # self.mean = float(-1) #test
        # self.var = float(-1) #test
        self.alpha = float(-1)
        self.bata = float(-1)
        
class Favre():
    def __init__(self):
        self.velo =  Velo()
        self.mixfrc = Mixfrc()
        
class Mean():
    def __init__(self):
        self.velo =  Velo()
        self.mixfrc = Mixfrc()
        
class Std():
    def __init__(self):
        self.velo =  Velo()
        self.mixfrc = Mixfrc()

class Var():
    def __init__(self):
        self.velo =  Velo()
        self.mixfrc = Mixfrc()
        
class Gridcell(): #细胞网格类型
    def __init__(self):          
        self.state='undefined' #分为四种 'undefined' 'wall' 'cold' 'burnt'分别表示初始化 不感兴趣 冷态 已经燃烧
        self.coord = Coord()
        self.velo = Velo()
        self.mixfrc = Mixfrc()
        self.favre = Favre()
        self.mean = Mean()
        self.std = Std()
        self.var = Var()
        self.L_turb = 1E-3 #湍流积分尺度，单位 m
        self.evapo_rate = 0 #单位 kg*s^-1*m^-3
        self.dens = -1.0 #单位 kg*m^-3
        self.visc = -1.0 #单位 kg*m^-3
        self.TKE = -1.0
    
class Particle(): #火焰例子类型
    def __init__(self):
        self.coord = Coord()
        self.velo =  Velo()
        self.mixfrc =  Mixfrc()
        self.TKE = init_float
        self.dissip = init_float
        self.omega = init_float
        self.velo_std = init_float
        self.Ka = init_float
        self.kinematic_visc = init_float
    
    def getcell(self):#返回[m,n,k]
        _m = math.floor((self.coord.x-x_coord_start) / grid_dist)
        _n = math.floor((self.coord.y-y_coord_start) / grid_dist)
        _k = math.floor((self.coord.z-z_coord_start) / grid_dist)
        return _m,_n,_k

print('0: set global variables...')
        
        
'step 1: load interpolated result'
grid_points_valid   = np.load('data/grid_points_valid.npy')
grid_position_valid = np.load('data/grid_position_valid.npy')
grid_V_ABS_valid    = np.load('data/grid_V_ABS_valid.npy')
grid_U_valid        = np.load('data/grid_U_valid.npy')
grid_V_valid        = np.load('data/grid_V_valid.npy')
grid_W_valid        = np.load('data/grid_W_valid.npy')
grid_CH4_valid      = np.load('data/grid_CH4_valid.npy')
grid_O2_valid       = np.load('data/grid_O2_valid.npy')
grid_DENS_valid     = np.load('data/grid_DENS_valid.npy')
grid_VISC_valid     = np.load('data/grid_VISC_valid.npy')
grid_TED_valid      = np.load('data/grid_TED_valid.npy') #Turbulence Eddy Dissipation
grid_TKE_valid      = np.load('data/grid_TKE_valid.npy') #Turbulence Kinetic Energy

grid_num = np.size(grid_points_valid,0)
print('1: cfx result loaded...')


''' varibles '''

loop_iters = 20
time_iters = 200

projectfolder = './Ka' + str(Ka_crit) +'_'
if (False == os.path.exists(projectfolder)) : 
    os.mkdir(projectfolder)


''' SET ZONE'''
ignit_x_min = 0e-3
ignit_x_max = 5e-3#55e-3
ignit_x_gap = 5e-3
ignit_x_num = math.floor((ignit_x_max-ignit_x_min)/ignit_x_gap) + 1

ignit_z_min = 0e-3
ignit_z_max = 5e-3#25e-3
ignit_z_gap = 5e-3
ignit_z_num = math.floor((ignit_z_max-ignit_z_min)/ignit_z_gap) + 1

probability_map = np.zeros([ignit_x_num,ignit_z_num]) #record ignition process

for ignit_x_iter in range(0,ignit_x_num):
    for ignit_z_iter in range(0,ignit_z_num):
       
        
        strtime = time.strftime('%m%d_%H%M')
        
        pi_ign_list = np.zeros([time_iters,loop_iters]) #record ignition process
        
        success_count = 0  #count the success events ignition position

        kernel_radius = 2E-3 
        kernel_x = ignit_x_min + ignit_x_iter * ignit_x_gap   #16    
        kernel_y = 0E-3     
        kernel_z = ignit_z_min + ignit_z_iter * ignit_z_gap  #7
        
        foldername = projectfolder+'/' + str(int(1000 * kernel_x)).zfill(2)+'mm_' + str(int(1000 * kernel_z)).zfill(2)+'mm_' #+ strtime
        if (False == os.path.exists(foldername)) : 
            os.mkdir(foldername)
        
        t_start = time.time()

        ''' Loop 1: MonteCarlo loop '''
        for loop_iter in range(loop_iters):
            
            savepath = foldername +'/'+ str(loop_iter).zfill(5)
            if (False == os.path.exists(savepath)) : 
                os.mkdir(savepath)

            total_count = 0 #最初冷态细胞个数
            burnt_count = 0 #目前已燃细胞个数
            burntlist = np.array([],dtype = float).reshape([0,6])
            pi_ign = float(0)
        
            
            'step 2: initialize a 3D grid contain the coordinates and states'
            x_num = math.ceil(x_size / grid_dist)
            y_num = math.ceil(y_size / grid_dist)
            z_num = math.ceil(z_size / grid_dist)
            
            grid = np.zeros((x_num, y_num, z_num),dtype=object) 
            
            count = 0     
            coldcount = 0 
            wallcount = 0
                            
            for m in range(x_num):
                for n in range(y_num):
                    for k in range(z_num):         
                        #initialize Gridcell type
                        grid[m][n][k]=Gridcell() #初始化类型为Gridcell
                        #set coordinates
                        grid[m][n][k].coord.x, grid[m][n][k].coord.y, grid[m][n][k].coord.z\
                        = m*grid_dist + x_coord_start, n*grid_dist + y_coord_start, k*grid_dist + z_coord_start 
                        grid[m][n][k].state = 'wall'
                        wallcount = wallcount + 1
                        
            for point in grid_position_valid:
                grid[point[0]][point[1]][point[2]].state = 'cold'
                wallcount = wallcount - 1
                coldcount = coldcount + 1
                        
            total_count = coldcount
            print('2: uninitialized:',x_num*y_num*z_num - coldcount - wallcount,' cells')
            print('   coldstate: ',coldcount,' cells')
            print('   walls: ',wallcount,' cells')
            
            
            
            # test = grid_U_valid[:,1]**2 +grid_V_valid[:,1]**2 +grid_W_valid[:,1]**2 
            # test1[:] = test[:]**0.5 - grid_V_ABS_valid[:,1]
            # np.std(test1)
            # np.std(grid_V_ABS_valid)
            
            
            'step 3：initialize mixture fraction / velocity / density / viscosity /turbulence kinetic energy'
            #get favre-averaged-velocity
            grid_FAVRECH4_valid = grid_CH4_valid[:,0] 
            grid_FAVREU_valid = grid_U_valid[:,0] 
            grid_FAVREV_valid = grid_V_valid[:,0] 
            grid_FAVREW_valid = grid_W_valid[:,0] 
            
            count = 0
            for i in range(grid_num):
                count = count + 1
                
                #get statistic value 
                m,n,k = grid_position_valid[i]
            
                grid[m][n][k].mean.mixfrc.value = grid_CH4_valid[i,0]
                grid[m][n][k].std.mixfrc.value = grid_CH4_valid[i,1]
                grid[m][n][k].var.mixfrc.value = grid_CH4_valid[i,1]**2 
            
                grid[m][n][k].dens = grid_DENS_valid[i]
                grid[m][n][k].visc = grid_VISC_valid[i]
                grid[m][n][k].TKE = grid_TKE_valid[i]  
                
                grid[m][n][k].mean.velo.x = grid_U_valid[i,0]
                grid[m][n][k].mean.velo.y = grid_V_valid[i,0]
                grid[m][n][k].mean.velo.z = grid_W_valid[i,0]
                
                # grid[m][n][k].std.velo.x = grid_U_valid[i,1]
                # grid[m][n][k].std.velo.y = grid_V_valid[i,1]
                # grid[m][n][k].std.velo.z = grid_W_valid[i,1]
                
                #各项同性假设
                grid[m][n][k].std.velo.x = ((2/3) *grid[m][n][k].TKE)**0.5 
                grid[m][n][k].std.velo.y = ((2/3) *grid[m][n][k].TKE)**0.5
                grid[m][n][k].std.velo.z = ((2/3) *grid[m][n][k].TKE)**0.5
                
                #test m,n,k = 75,55,10 m,n,k = 55,55,0
                
                
                #generate mixture fraction using random "BETA DISTRIBUTION" if no error
                if(grid[m][n][k].var.mixfrc.value != 0 and grid[m][n][k].mean.mixfrc.value != 0): 
            
                    grid[m][n][k].mixfrc.alpha = (1 - grid[m][n][k].mean.mixfrc.value) * grid[m][n][k].mean.mixfrc.value**2 / grid[m][n][k].var.mixfrc.value - grid[m][n][k].mean.mixfrc.value  
                    grid[m][n][k].mixfrc.beta = grid[m][n][k].mixfrc.alpha*(1/grid[m][n][k].mean.mixfrc.value - 1)
                    if(grid[m][n][k].mixfrc.alpha>0 and grid[m][n][k].mixfrc.beta>0):
                        grid[m][n][k].mixfrc.value = np.random.beta(grid[m][n][k].mixfrc.alpha, grid[m][n][k].mixfrc.beta)
                    else:
                        grid[m][n][k].mixfrc.value = grid[m][n][k].mean.mixfrc.value
                else:
                    grid[m][n][k].mixfrc.value = grid[m][n][k].mean.mixfrc.value
                    
                
                #apply favre-averaged velocity & favre-averaged mixture fraction
                grid[m][n][k].favre.mixfrc.value = grid_FAVRECH4_valid[i]
                grid[m][n][k].favre.velo.x = grid_FAVREU_valid[i]
                grid[m][n][k].favre.velo.y = grid_FAVREV_valid[i]
                grid[m][n][k].favre.velo.z = grid_FAVREW_valid[i]
                
                #generate velocity using random "GAUSS DISTRIBUTION"
                grid[m][n][k].velo.x = random.gauss(grid[m][n][k].mean.velo.x , grid[m][n][k].std.velo.x) #网格后续只有favre averaged velocity参与
                grid[m][n][k].velo.y = random.gauss(grid[m][n][k].mean.velo.y , grid[m][n][k].std.velo.y)
                grid[m][n][k].velo.z = random.gauss(grid[m][n][k].mean.velo.z , grid[m][n][k].std.velo.z)
            
                            
            print('3: mixfrction and velocity initialized:',count,'cells')
            # print('\n') 
            
            # grid[7][6][3].mixfrc.value
            # grid[7][6][4].mixfrc.mean
            # grid[6][6][5].mixfrc.var
            
            
            
            'step 4: initialize the flame kernel and the flame particles'   
            kernel_count = 0
            flame_particles=list()#用来存放火焰粒子
            
            for m in range(x_num):
                for n in range(y_num):
                    for k in range(z_num):
                        check = ((grid[m][n][k].coord.x-kernel_x)**2 + (grid[m][n][k].coord.y-kernel_y)**2 + (grid[m][n][k].coord.z - kernel_z)**2) < kernel_radius**2
                        if (check == True):
                            if('cold' == grid[m][n][k].state):
                                grid[m][n][k].state = 'burnt'
                                burnt_count += 1
                                burntlist = np.append(burntlist,np.array([[m,n,k,grid[m][n][k].coord.x,grid[m][n][k].coord.y,grid[m][n][k].coord.z]]),axis=0)
                                
                                kernel_count += 1
                                new = Particle()
                                new.coord = copy.deepcopy(grid[m][n][k].coord)
                                new.velo = copy.deepcopy(grid[m][n][k].velo)
                                new.mixfrc = copy.deepcopy(grid[m][n][k].mixfrc)
                                new.ka= -1
                                new.velo_std = np.linalg.norm([grid[m][n][k].std.velo.x,grid[m][n][k].std.velo.y,grid[m][n][k].std.velo.z])
                                new.TKE = copy.deepcopy(grid[m][n][k].TKE)
                                
                                flame_particles.append(new)
                     
            print('4: flame kernel initialized:',kernel_count,'cells')
            pi_ign = float(burnt_count) / float(total_count)   
            print('   PI_ign = ',pi_ign*100,'%')            
            
            
            
            
            ''' Loop 2: timestep update loop '''
            for timestep in range(time_iters):
                
                'step 5: update the position of flame particles at t=t+dt '
                '''
                更新了位置，并剔除了不在拷问区域的粒子
                '''
                i=0 #表示判断次数
                j=0 #表示目前flame_particles的第j个
                current_amount = len(flame_particles)
                while(i<current_amount):
                    #update coordinates
                    flame_particles[j].coord.x += flame_particles[j].velo.x * dt #单位 m = m/s * s
                    flame_particles[j].coord.y += flame_particles[j].velo.y * dt #单位 m = m/s * s
                    flame_particles[j].coord.z += flame_particles[j].velo.z * dt #单位 m = m/s * s
                    
                    #if particle went in the wall, delete them
                    check = False #不符合后续条件则False，待删除
                    m,n,k = flame_particles[j].getcell()
                    #在程序几何范围内
                    if(0<=m<x_num and 0<=n<y_num and 0<=k<z_num):
                        #‘非wall’部分则Ture
                        if(grid[m][n][k].state !='wall'):
                            check = True
                    #不在几何范围或者进入'wall'则删除
                    if(check == False):
                        del flame_particles[j]
                        j -= 1
                    #更新i和j   
                    i += 1
                    j += 1
                
                # current_amount = len(flame_particles)
                # if(not current_amount):
                #     print('没有火焰粒子\n')
                # else:
                #     print('打印第0个粒子坐标(test)')
                #     print(flame_particles[0].coord.x,flame_particles[0].coord.y,flame_particles[0].coord.z)
                
                # print('5: update the position of flame particles at t=t+dt')
                
                
                
                
                'step 6: update the velocity and mixfrction of flame particles at t=t+dt'
                i=0
                j=0
                current_amount = len(flame_particles)
                while(i<current_amount):
                    #get current cell
                    m,n,k = flame_particles[j].getcell()
                    
                    #calculate velostd 计算湍流脉动速度
                    # flame_particles[j].velo_std = np.linalg.norm([grid[m][n][k].std.velo.x,grid[m][n][k].std.velo.y,grid[m][n][k].std.velo.z])
                    flame_particles[j].velo_std = ((2/3)*grid[m][n][k].TKE) **0.5 
                    
                    #calculate turbulence kinetic energy 计算湍流动能
                    # flame_particles[j].TKE = (3/2) * ((flame_particles[j].velo_std)**2) #单位 (m^2/s^2)
                    flame_particles[j].TKE = grid[m][n][k].TKE
                    
                    #calculate inverse integral turblance timescale 计算反向积分湍流时间尺度
                    flame_particles[j].omega = flame_particles[j].velo_std / grid[m][n][k].L_turb   #单位 (m/s) / (m) = (1/s)
                    
                    #calculate turbulence dissipation 计算湍流动能耗散率
                    flame_particles[j].dissip = flame_particles[j].TKE * flame_particles[j].omega #单位(m^2/s^2) * (1/s) = (m^2/s^3)
                    
                    #update flame particle velocity 更新火焰粒子速度
                    flame_particles[j].velo.x += \
                        dt * ( -(0.5+0.75*C0)*flame_particles[j].omega * (flame_particles[j].velo.x - grid[m][n][k].favre.velo.x ) )\
                        + (C0 * flame_particles[j].dissip * dt)**0.5 * random.gauss(0,1)
                    flame_particles[j].velo.y += \
                        dt * ( -(0.5+0.75*C0)*flame_particles[j].omega * (flame_particles[j].velo.y - grid[m][n][k].favre.velo.y ) )\
                        + (C0 * flame_particles[j].dissip * dt)**0.5 * random.gauss(0,1)
                    flame_particles[j].velo.z += \
                        dt * ( -(0.5+0.75*C0)*flame_particles[j].omega * (flame_particles[j].velo.z - grid[m][n][k].favre.velo.z ) )\
                        + (C0 * flame_particles[j].dissip * dt)**0.5 * random.gauss(0,1)
                    
                    
                    #update flame particle velocity 更新火焰粒子混合分数
                    flame_particles[j].mixfrc.value += (-0.5)* C1 * flame_particles[j].omega * (flame_particles[j].mixfrc.value - grid[m][n][k].favre.mixfrc.value ) * dt \
                        + (1-flame_particles[j].mixfrc.value) * grid[m][n][k].evapo_rate * dt / grid[m][n][k].dens
                    
                    i += 1
                    j += 1
                    
                # current_amount = len(flame_particles)
                # if(not current_amount):
                #     print('没有火焰粒子\n')
                # else:
                #     j=2
                #     print('打印第',j,'个粒子坐标(test)')
                #     print('湍流脉动速度: ',flame_particles[j].velo_std,'(m/s) \n',\
                #           '湍流动能: ',flame_particles[j].kinetic,'(m^2/s^2) \n', \
                #           '反向积分湍流时间尺度: ',flame_particles[j].omega,'(1/s) \n', \
                #           '湍流动能耗散率: ',flame_particles[j].dissip,'(m^2/s^3) \n',\
                #           '新的速度\n','x:',flame_particles[j].velo.x,'\ny:',flame_particles[j].velo.y,'\nz:',flame_particles[j].velo.z,'\n'\
                #           '新的混合分数：',flame_particles[j].mixfrc.value,'\n'
                #           )
                # if(not current_amount):
                #     print('没有火焰粒子\n')
                # else:
                #     j=4
                #     print('打印第',j,'个粒子坐标(test)')
                #     print('湍流脉动速度: ',flame_particles[j].velo_std,'(m/s) \n',\
                #           '湍流动能: ',flame_particles[j].kinetic,'(m^2/s^2) \n', \
                #           '反向积分湍流时间尺度: ',flame_particles[j].omega,'(1/s) \n', \
                #           '湍流动能耗散率: ',flame_particles[j].dissip,'(m^2/s^3) \n',\
                #           '新的速度\n','x:',flame_particles[j].velo.x,'\ny:',flame_particles[j].velo.y,'\nz:',flame_particles[j].velo.z,'\n'\
                #           '新的混合分数：',flame_particles[j].mixfrc.value,'\n'
                #           )
                
                
                # print('6: update of the velocity and mixfrction of flame particles')
                
                
                
                
                'step 7: generate new flame particles'
                i=0
                j=0
                current_amount = len(flame_particles)
                while(i<current_amount):
                    #get current cell
                    m,n,k = flame_particles[j].getcell()
                    
                    check = (grid[m][n][k].state == 'cold')
                    if (check):
                        grid[m][n][k].state = 'burnt'
                        burnt_count += 1
                        burntlist = np.append(burntlist,np.array([[m,n,k,grid[m][n][k].coord.x,grid[m][n][k].coord.y,grid[m][n][k].coord.z]]),axis=0)
                        
                        new = Particle()
                        new.coord = copy.deepcopy(grid[m][n][k].coord)
                        new.velo = copy.deepcopy(grid[m][n][k].velo)
                        new.mixfrc = copy.deepcopy(grid[m][n][k].mixfrc)
                        new.ka= -1
                        new.TKE = copy.deepcopy(grid[m][n][k].TKE)
                        new.velo_std = np.linalg.norm([grid[m][n][k].std.velo.x,grid[m][n][k].std.velo.y,grid[m][n][k].std.velo.z])
                        flame_particles.append(new)
                        
                    i += 1
                    j += 1
                
                # print('7: generate new flame particles')
                
                
                
                
                
                'step 8: kill flame particles'
                # """
                # CH4 按照文献获得混合分数对应的层流火焰速度 + 运动粘度（粘度/密度） + 湍流脉动速度 ==> Ka 
                # """
                # stoi = 16 / (64/0.232) #stoichiometric 化学计量燃料和空气质量比
                # x = np.linspace(0.6, 1.5, num=10)
                # y = np.array([8,16,24.5,30,35,37,34,23,14,10],dtype=float) #单位 cm/s
                # y = y * 1E-2 #单位 m/s
                # coeff = np.polyfit(x, y, 4)
                # func_LFS = np.poly1d(coeff)
                # # xtest = np.linspace(0.6, 1.5, num=100)
                # # print(coeff)
                # # p = plt.plot(x, y, 'rx')
                # # p = plt.plot(xtest, f(xtest))
                eq_stoi = 44 / (160/0.232) #stoichiometric 化学计量燃料和空气质量比
        
                
                i=0
                j=0
                current_amount = len(flame_particles)
                test = np.zeros([current_amount])
                while(i<current_amount):
                    #get current cell
                    m,n,k = flame_particles[j].getcell()
                    
                    kinematic_visc = grid[m][n][k].visc / grid[m][n][k].dens
                    
                    eq_ratio = 1 / ((1/flame_particles[j].mixfrc.value - 1) * eq_stoi)  
                    
                    LFS = (1E-2) * (34.22 + -138.65 * (eq_ratio - 1.08)**2) #mnk 73 48 10 .coord.x grid[m][n][k].coord.y grid[m][n][k].coord.z
                    if(LFS<=1E-6):
                        LFS = 1E-6
        
                    flame_particles[j].Ka = 0.157 * (kinematic_visc * flame_particles[j].velo_std**3 / grid[m][n][k].L_turb)**0.5 / LFS
                    test[i] = flame_particles[j].Ka
                    #flame_particles[50].getcell()  # m,n,k=79, 55, 7
                    
                    if(flame_particles[j].Ka > Ka_crit ):
                        del flame_particles[j]
                        j=j-1
                    
                    i += 1
                    j += 1
                
                
                
                # print('8: claculate quenching criterion and kill old flame particles')
                
                
                
                
                
                'step 9: calculate the current ignition progress Pi_ign'
                current_amount = len(flame_particles)
                # print('   fparticle_amount = ',current_amount) 
                # if(current_amount == 0):
                    # print('The end of a progress')
                    # print('  ignition time:',)
                    
                pi_ign = float(burnt_count) / float(total_count)   
                print('  Loop iter:' + str(loop_iter) + '  Time step:' + str(timestep) + '  PI_ign = ',pi_ign*100,'%') 
                pi_ign_list[timestep][loop_iter] = pi_ign
                
            
            
                'for every 5 step'
                if(timestep%5 == 0):
                    fname = str(timestep)
                    fname = fname.zfill(5)
                    
                    'print the burnt cell'
                    # fig = plt.figure() # 创建一个画布figure，然后在这个画布上加各种元素。
                    # ax = Axes3D(fig)
                    # ax.scatter(burntlist[:,3],burntlist[:,4],burntlist[:,5]) # 画出(xs1,ys1,zs1)的散点图。
                    # ax.set_xlabel('X') # 画出坐标轴
                    # ax.set_ylabel('Y')
                    # ax.set_zlabel('Z')
                    # ax.set_xlim(x_coord_start,x_coord_start+x_size )
                    # ax.set_ylim(y_coord_start,y_coord_start+y_size )
                    # ax.set_zlim(z_coord_start,z_coord_start+z_size )
                    # ax.azim = -143
                    # ax.elev = 0  
                    # # plt.show()
                    # # plt.savefig('./pics'+foldertime+'/'+fname+'.png')
                    # plt.savefig(savepath+'/'+fname+'.png')
                    # plt.close(fig)
                
                    'save the burnt cell axis to csv'
                    colname = ['m','n','k','x','y','z']
                    cell_dataframe = pd.DataFrame(data = burntlist,columns = colname)
                    cell_dataframe.to_csv(savepath+'/'+fname+'.csv',index=None)
                    
            
            #if 
            if(pi_ign*100 > Pi_crit): success_count = success_count + 1
            
        'possibility'    
        possibility = success_count / loop_iters
        probability_map[ignit_x_iter,ignit_z_iter] = possibility
        
        
        'Pi'
        np.save(foldername+'/'+'pi_ign_list',pi_ign_list)
        pi_listx = np.arange(time_iters)
        
        fig = plt.figure()
        for i in range(loop_iters):  
            pi_listy = pi_ign_list[:,i]*100
        
            plt.plot(pi_listx, pi_listy, label = 'no: %d'%i)
            plt.legend(loc='upper left')
            
        # ax.set_xlabel('X') # 画出坐标轴 
        # ax.set_ylabel('Y') 
        plt.axis([-1,time_iters,0,1.2*np.max(100*pi_ign_list)])
        plt.xlabel('timestep, delta = 5E-3 (s)')
        plt.ylabel('burnt cell %')
        plt.title('possibility = %.2f%% when Ka = %.2f' %(100*possibility,Ka_crit))
        
        # plt.show()
        
        plt.savefig(foldername +'/'+'PI.png')
        
        
        'timer'
        t_end = time.time()
        print('loop iters:%d'%loop_iters)
        print('time iters:%d'%time_iters)
        print('ellapsed time:','%dm'%math.floor((int(t_end)-int(t_start))/60),'%ds'%int((t_end-t_start)%60))


        
possibility_df = pd.DataFrame(data = probability_map,columns=None)
possibility_df.to_csv(projectfolder+'/'+'P_map.csv',index=None,header = False)


















            